/***********************************************
*INVERSE COVARIANCE WEIGHTED INDEX
*This ado file creates an index using inverse covariance weighting, a la Anderson (2008).
*It was written by Cyrus Samii, NYU (cds2083@nyu.edu)
*For an updated version, see http://cyrussamii.com/?p=2656

*The command is 'make_index' and it takes three arguments in a required sequence 
*	(1) The prefix for the name of the new index
*	(2) The sampling weight (used in the covariance matrix)
*	(3) The name of the local macro with the variables for the weighted index.
	
*For example: 
*	make_index recid_direct `wgt' `recid_direct' 
	
*produces the index recid_direct_index using the weight in `wgt' and the index vars in local
*macro `recid_direct'

*References
*
*Anderson ML (2008) "Multiple Inference & Gender Differences..." JASA 103(484).

***********************************************/

capture program drop make_index
program make_index
version 11.1
    syntax anything [if]
    gettoken newname anything: anything
    gettoken wgt anything: anything
    local Xvars `anything'
	marksample touse
	quietly: cor `Xvars' [aweight=`wgt']
	matrix S = r(C)
  	mata: icwxmata(("`Xvars'"),"`wgt'", "index")
	rename index index_`newname'
end
 
mata:
	mata clear
	function icwxmata(xvars, wgts, indexname)
	{
		st_view(X=0,.,tokens(xvars))
		st_view(wgt=.,.,wgts)
		nr = rows(X)
		nc = cols(X)
		wgtst = wgt/sum(wgt)
		wgtstdM = J(1,nc,1) # wgtst
		wgtdmeans = colsum(X:*wgtstdM)
		wgtdmeandevs = (X - J(nr,1,1) # wgtdmeans)
		wgtdstds = sqrt(colsum(wgt:*(wgtdmeandevs:*wgtdmeandevs)):/(sum(wgt)-1))
		Xs = wgtdmeandevs:/(J(nr,1,1) # wgtdstds)
		invS = invsym(st_matrix("S"))
		ivec = J(nc,1,1)
		indexout_sc = (invsym(ivec'*invS*ivec)*ivec'*invS*Xs')'
		indexout = (indexout_sc - J(nr,1,1)# mean(indexout_sc))/sqrt(variance(indexout_sc))
		st_addvar("float",indexname)
		st_store(.,indexname,indexout)
	}
end
